Non-equilibrum dynamics in the strongly excited inhomogeneous Dicke model 
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Using the exact eigenstates of the inhomogeneous Dicke model obtained by numerically solving 
the Bethe equations, we study the decay of bosonic excitations due to the coupling of the mode to 
an ensemble of two-level (spin 1/2) systems. We compare the quantum time-evolution of the bosonic 
mode population with the mean-field description which, for a few bosons agree up to a relatively 
long Ehrenfest time. We demonstrate that additional excitations lead to a dramatic shortening of 
the period of validity of the mean-field analysis. However, even in the limit where the number of 
bosons equal the number of spins, the initial instability remains adequately described by the mean- 
field approach leading to a finite, albeit short, Ehrenfest time. Through finite size analysis, we 
also present indications that the mean-field approach could still provide an adequate description for 
thermodynamically large systems even at long times. However, for mesoscopic systems one cannot 
expect it to capture the behavior beyond the initial decay stage in the limit of an extremely large 
number of excitations. 

PACS numbers: 78.20.Bh,42.50.Pq,02.30.Ik 



I. INTRODUCTION 

Since they constitute one of the broad classes of 
proposed physical realizations of quantum computing 
devices 1-4 , coherently interacting light-matter systems 
have received lately a considerable level of attention. The 
interest has been further enhanced by relatively recent 
progress in a variety of systems ranging from polaritons 
in quantum wells 5-7 to semiconductor quantum dots 8-10 
which nowadays make it possible to realize solid-state 
based quantum systems which couple to a single photon 
eigenmode of optical microcavities. At the same time the 
many-body effects in these systems are becoming increas- 
ingly important in the context of engineering of semicon- 
ductor lasers 11-14 motivated by potentially a great en- 
hancement in performance of designs in which quantum 
dots serve as an active medium 15 . 

While an ideal system of identical two- level (spin 1/2) 
emitters coupled uniformly to a single light-mode is de- 
scribable in terms of the Dicke Hamiltonian 16,17 , a more 
realistic setup would need to include possible inhomo- 
geneities in both the energy splitting of the individual 
spins and in their respective coupling strengths to the 
bosonic light mode. In previous studies 19 ' 20 , the com- 
parison of the quantum and mean-field dynamics of the 
resulting generalized Dicke model 

N N 

j'=i i=i 

was performed by solving the Schrodinger equation in the 
limit of small excitation numbers. Since the total number 
of excitations M — tfb + J2i (Sf + \) is conserved, the 
dimension of the Hilbert space involved in the unitary 



evolution of the system is drastically reduced. Thus it 
became possible to solve the explicit time dependence of 
every quantum amplitude involved. 

In this work, we revisit this problem by exploiting 
the quantum integrability of a certain generalized Dicke 
model. Using the algebraic Bethe Ansatz (ABA), one 
can numerically compute exact eigenstates of the system 
and study its dynamics rendered trivial by the use of the 
proper eigenbasis. Additionally, this approach allows, in 
the strong coupling regime, a drastic truncation of the 
Hilbert space granting access to relatively large system 
sizes. On the other hand, integrability imposes the con- 
straint that every spin-like subsystem be uniformly cou- 
pled to the bosonic mode and consequently, we study the 
specific Hamiltonian 

JV N 

J'=l i=i 

It is known 19 that the eigenstates of Eq. (1) are the 
eigenstates of Eq. (2) with V = ^T,f=i V?/N, at least 21 
when the number of excitations is small M -C N. 

As in Ref. [20], this work focuses on the decay of 
a number of bosonic excitations due to the coupling of 
the mode to an initially unexcited set of two-level emit- 
ters. Such a state could, in principle, be obtained by first 
preparing the spin system in its ground-state then excit- 
ing the bosonic mode via an external radiation pulse. For 
small number of excitations, a crossover between two dis- 
tinct regimes was found. At weak coupling, the bosonic 
occupation number (b^b) undergoes an exponential decay 
and at strong enough coupling, the occupation undergoes 
non-decaying periodic oscillation with frequency which is 
enhanced by the Dicke supperradiance effect. This is a 
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dynamical counterpart of the Hepp and Lieb superradi- 
ance quantum phase transition. 23 

Starting from a large number of excitations, for weak 
coupling strength, the dimension of the necessary Hilbert 
space severely limits the system sizes treatable using the 
ABA. Still, for very small systems, we do find rapidly 
decaying short-time dynamics as in the mean-field treat- 
ment. However, our capacity to make quantitative com- 
parisons with the mean-field approach is hindered since 
its validity is necessarily limited to large systems. Con- 
sequently, the bulk of our results focuses on the strong 
coupling regime, where heavy truncation of the Hilbert 
space allows a nearly exact treatment of larger systems. 
In this case, for a small number of initial excitations, the 
spectrum obtained in the full quantum treatment is char- 
acterized by set of equally spaced frequencies, leading to 
nearly periodic real-time dynamics. The mean-field ap- 
proach, which also leads to periodical oscillating bosonic 
populations, then reproduces the quantum dynamics up 
to some relatively long Ehrenfest time at which both so- 
lutions start to differ significantly. 

When the number of initial excitations becomes of the 
order of the system size, the spectrum shows strong devi- 
ations from a harmonic progression no longer reproducing 
the periodic oscillations obtained in mean-field. Nonethe- 
less, we find that the mean-field approach remains valid 
for the very short-time dynamics of the system leading 
to a finite, but considerably shortened, Ehrenfest time. 

The crossover to the regime of periodic oscillations oc- 
curs when the superradiantly enhanced coupling with the 
boson mode becomes larger than the bandwidth of the 
spin energy splittings. For example, this effect can man- 
ifest itself as a suppression of the inhomogeneous broad- 
ening in a system of self-assembled quantum dots in a 
zero dimensional cavity. A realistic high in-plane density 
of InGaAs dots 24 ~ 5- 10 10 cm~ 2 will enhance coupling of 
a single dot to a photonic crystal cavity, 9 a micropillar, 25 
or a concave microcavity, 26 which is (<~ 3 — 100 fieV), by 
two to three orders of magnitude in the optical domain 27 
that is well inside or even above the range of natural 
bandwidths of ensembles of such dots 28 (^5 — 50 meV) . 
Thus for sufficiently high densities a spectrally broad- 
ened ensemble of the self-assembled dots will exchange 
all of its excitations with a boson mode without a decay, 
at least on a finite time scale, in the same way as an 
ideal atomic system without any broadening, making it 
suitable to engineer a high power semiconductor laser. 

When such a system is in the strong coupling regime, 
the effect of quantum fluctuations could be observed in a 
direct time-resolved measurement at different excitation 
powers. At low to intermediate number of excitations the 
many Rabi oscillations would be visible without any sig- 
nificant decay. However, when the number of excitations 
reaches <~ 80% of the number of spins, a strong decay 
on a time scale of a single period would appear due to 
quantum fluctuations. The dependence of this effect on 
the number of spins can be used to discriminate it from 
other sources of relaxation that occur at large probing 



powers such as charge or phonon fluctuations. 

The paper is organized as follows. Section II describes 
the exact solution (ABA) and the numerical techniques 
used to exploit it in order to compute the non-equilibrium 
dynamics we are interested in. Section III covers the 
mean-field approach to the same problem. The resulting 
behavior which stems from both descriptions are then 
compared and analyzed in Section IV where both the 
spectrums and the real-time dynamics are studied. 

II. THE EXACT SOLUTION 

In the following analysis, we use ej (single spin excita- 
tion energies) which are uniformly distributed e,_|_i — ej — 
ejN/ (N — 1) within a band of total width A = e^v — £i = 
N, here e<j is a "level spacing" for spins. This width A 
will serve as a natural energy scale of the problem. Fur- 
thermore, we introduce the Rabi frequency CI = V\/N, 
which is the oscillation frequency for bosons in the case 
of equal splittings = e. The ratio Cl/A is thus a di- 
mensionless parameter that we will use to specify the 
coupling strength. 

A. Constructing Eigenstates 

We exactly solve the Dicke Hamiltonian (2) using the 
method introduced in Rcfs. [30] and [31]. Hereby, we ex- 
ploit the quantum integrability of the Dicke Hamiltonian, 
which was proven in [37]. Unnormalized eigenstates can 
thus be constructed by creating M pseudo-particles, 

M 

nS + (A Q )|0;|...|), (3) 

a=l 

on the vacuum state |0;i . . . \.) which contains no bosons 
and has all spins in their lowest energy states. For the 
Dicke model (2) the creation operator takes the form 

N V 

S+(A Q ) = 6t + ^-— -S+. (4) 

J=l 

States of the form (3) become eigenstates of the Dicke 
Hamiltonian when the M rapidities A a fulfill the M Bethc 
equations 

M T/ \ i ^ t T 

\ - V U Aa ly. V 

which one can obtain from a straightforward application 
of the Hamiltonian (2) to a general state (3) with un- 
specified rapidities A a 32 . 

Being completely defined by a set of M rapidities {A}, 
we denote the (unnormalized) eigenstates by |{A}). The 
corresponding eigenenergies are simply given by 

M 

£{A} = E A - ( fi ) 

a=l 
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Due to numerical instabilities when trying to solve the 
Bethe ansatz Eqs. (5) directly (see for example [29]), we 
use the change of variables proposed in [30,31]. First we 
introduce the complex polynomial 



with 



M 



m = E 



a=l 



A 



(7) 



which we evaluate at the Zeeman splittings z = ej to 
obtain N new variables Kj = V 2 T(ej). For a set of 
rapidities that are a solution to (5) one can show that 
the corresponding new variables obey a set of quadratic 
equations 

v 2J2^k + V 2 M = K j (e j -u) + K*. (8) 
t—i j 

Note that instead of M, we now have N equations. As 
long as M < N, these equations are equivalent to the 
Bethe ansatz Eqs. (5), but they lack the numerical prob- 
lems mentioned before. For every quasi-particle number 
M, which is conserved in time, the equations allow for 
several solutions {K}, all of them in one to one corre- 
spondence to a given set of rapidities {A} and thus to a 
single cigenstate of the system. 

At V = 0, the eigenstates of the system are obviously 
Fock states with a definite number of spin excitations 
and bosons. For example, for N = 2 and M = 2 the 
Fock states are [2; 44), |l;t4) 5 |l;4-t) and |0;tt) where, 
in this notation, the number counts the bosonic excita- 
tions whereas the arrows represent the spin states. From 
the Bethe Eqs. (5) it immediately follows that the set of 
rapidities {A} of a Fock state consists of one A Q = uo for 
every boson and one A a = if the ith spin is excited. 
This leads to Ki = u> — if the ith spin is excited and 
Ki = if it is not. The total number of excitations M 
then differentiates between states with identical spin con- 
figurations but different number of bosons. For example, 
the state |l;t4) has {A} = {w,ei} and {K} = {0,w-ei}. 

For a desired final coupling V — Vj, we obtain the 
solutions {K} by deforming the solutions at V — by 
a stepwise increasing of V. As detailed in Rcf. [30], 

one can compute easily the nth first derivatives d ■ 
This provides an good initial guess for the solution at 
V + dV using the values Ki(V) through the truncated 
Taylor expansion. One can the refine this guess using 
a simple iterative Newton-Raphson algorithm applied to 
the quadratic system of Eqs. (8). The process can then 
be repeated until the target coupling value Vf is reached. 

Since the rapidities themselves are used to calculate 
physical quantities, we need to extract the set of M values 
{A} corresponding to a given set {K}. In this work it is 
achieved by using the fact that 



r(z) = 



Q(z) 
Q'(zV 



(9) 



M 



M 



Q(z)=l[(z-\ fj ) = J2Q a z 



M—a 



a=0 



The coefficients of this polynomial 



M 



Qa — (—1)" E ^fei " ' ^k a 



(10) 



(11) 



ki=l 



are the elementary symmetric polynomials which can be 
found by solving the linear system: 



M 

E 



(M - a)e 



M-a-l _ f£i M-a 
3 y2 j 



q _ ^± e M_ M M-\ 



(12) 

The last task is to find the rapidities {A} as the roots 
of the polynomial Q(z) by its coefficients Q a , which one 
can do using a number of root finding algorithms. Using 
this method, we can, in principle, calculate the rapidities 
{A} characterizing every single eigenstate of the system. 



B. Obtaining the Bosonic Occupation Number 

Eventually, we are interested in the time dependent 
bosonic occupation number (tfb)(t) = {^)(i)\b^b\ip(t)) 
with initial state |Af;4 • • • 4)- Expanding in the normal- 
ized eigenbasis \4>i) = |{A}j)/-\/({A}j|{A}j}, we obtain 

d 

w)\tfbw)) = E< M; +---+ i&x&i 6t *#i> 

(^|M;4...4) e*^"^)*, (13) 

where we denote by d the dimension of the Hilbert space. 
The matrix elements occurring in (13), as well as the 
norm of the eigenstates ({A}i|{A}i) can be computed us- 
ing Slavnov's formula [33]. Provided the set {/i} fulfills 
the Bethe ansatz Eqs. (5), one can write its overlap with 
a generic state of the form (3) as the determinant of an 
M by M matrix 



(MI{A}> - 



with 

G n 



IL>z( A fc - A;)rifc<z(Mfc - W) 



A ^ + Ei— -+ 



detG, (14) 



V ^7 - A 



A/3 - H/3 

{\p -M«) 2 ' 



(15) 



For the norms of eigenstates we hence obtain, in the limit 
M -+ {A}, 



({A}|{A})=det(W), 



(16) 
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where 



2V 2 



(A/3 - X a ) 2 



(17) 



In a way similar to [34], we can furthermore derive a 
single determinant expression for the form factor appear- 
ing in Eq. (13): 



m\tfb\{\}) = 



n^( A i-M/=) 



IL>i( A fc - x i)U k <i(p-k - w) 

with the M x M matrices Q defined by 

(A/3 - A ; ) 



det(G + Q), (18) 



n 



(A/3 - w) ' 



(19) 



On the other hand, since \M;\. ... \) = 
l/VM! 6^ M |0;| ... |), we can rewrite the overlaps be- 
tween eigenstates and the initial state by noting that only 
the bosonic parts of the eigenstates (4) do not vanish in 
this projection 



{MU...i\4>i) = 



Ml 



<{A}|{A})' 



(20) 



Hence, we are left with the norms of the eigenstates 
({A}|{A}) which are computed by (17). In this fashion, 
we are able to compute the bosonic occupation (13) fully 
in terms of the rapidities of all eigenstates. 



C. Hilbert Space Truncation 

Although every term in Eq. (13) can be easily com- 
puted, the double sum over the complete Hilbert space 
remains remarkably large. For a system with O(10) spins 
this sum already becomes impossible to perform fully. 
In this work we rely on the fact that, at weak enough 
(CI <C A /AT) or strong enough (O 3> A) coupling, the 
main contributions to (13) comes from only a small num- 
ber of eigenstates. 

The truncation scheme for very small coupling works 
as follows. For V = 0, the initial state \M;l ... I) is 
an eigenstate itself and therefore the only relevant state 
for the given scenario. Perturbation theory then pro- 
vides a natural hierarchy such that states where a sin- 
gle excitation is swapped from the bosonic mode to a 
spin, are the most relevant ones. Therefore, keeping only 
states with a single spin excitation, for example the state 
\M — 1; t-i • • • 4-)> would lead to a large contribution. One 
could then add two spin-excitations states and so on. 

At strong coupling, considering Eqs. (3) and (4) we 
see that, as V — > oo, any finite rapidity leads to an ex- 
citation (created by S + (A)) which exclusively affects the 
spin sector. On the other hand, any rapidity, which di- 
verges when V — > oo, creates an excitation that signif- 
icantly populates the bosonic mode. When looking at 



the projection of any eigenstate onto the purely bosonic 
initial state Eq. (20), we can infer that, at large enough 
couplings, only the eigenstates for which every one of 
the M rapidities diverge will lead to significant overlaps. 
Since the form factors (cf>i\tfb\(f)j) in (13) are bounded 
by M, only the eigenstates with all rapidities diverging 
are needed in the V — > oo limit. At finite but large V, 
similar argument can be used to show that states with 
M — 1 diverging rapidities are the first ones to become 
important, and so on. This provides a rough ordering of 
the relative contribution to the sum of different classes 
of eigenstates. 

Having identified the heavily contributing states, we 
can truncate the sum in Eq. (13) to a smaller dimension 
d using only a subset of the most important eigenstates. 
We estimate the validity of the truncated sum by com- 
puting the sum rule 



E(d) = J2\(MU...i\4>i)f 



(21) 



i=l 



For the complete sum d = d we have T,(d) = 1 since we 
are projecting the initial states on the complete eigen- 
basis. In view of bounding the total error 5 in the 
time evolution for arbitrary times, we require at least 
S T (J) >1-S. While the factor b^) in Eq. (13) is 
different for different eigenstates, it is of the same order 
of magnitude for eigenstates containing the same num- 
ber of divergent rapidities, which can be seen from Eq. 
(4). Therefore, the saturation of the sum rule is a clear 
indication of the error generated by the truncation. Nev- 
ertheless, we crosscheck the validity of the truncated sum 
by checking if (tfb)(0)/M > 1 - 6 at t = 0. Due to the 
trivial time evolution in the true eigenbasis, the absolute 
error should remain bounded by this initial value. 



III. MEAN-FIELD ANALYSIS 

In this section we study the dynamics of the model Eq. 
(2) using the mean-field approximation. We derive the 
Hamilton equations of motion in terms of the expectation 
values of the boson and spin operators, and solve them 
for the initial conditions from Sec. II B for arbitrary M. 

Heisenberg equations of motion for quantum operators 
are obtained from Eq. (2) by use of commutation rela- 
tions, e.g. b = i [H, b]. The complete set of equations of 
evolution for the boson and spin operators is 



Sj — Bj x Sj 



N 



3 = 1 



(22) 
(23) 



where z and in-plane components of vector Bj = 
(2Vb x , 2Vb y , €j) are a single spin splitting and boson op- 
erators; b x + ib y — and b x — ib y — b. When the 
harmonic oscillator is highly excited the boson opera- 
tor can be approximated by a time-dependent c-number 
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(b) = a which makes the systems of Eqs. (22, 23) linear 
in operators. Here (...) is the time-dependent quantum 
mechanical expectation value. Averaging the linearized 
equations over an initial state we obtained the dynamical 
mean-field equations, 



— x Cj 

N 



(24) 
(25) 



where Cj = (Sj) is a set of N vectors of length \Cj\ = 1/2 
and Bj = (2Va x ,2Va y ,€j); 



here Cj = C* 



iCj and 



a = a x — ia y . 

Alternatively, Eqs. (24) and (25) can be derived us- 
ing Dirac's analogy for dynamical variables: Commuta- 
tion relations between quantum operators correspond to 
Poisson brackets between classical degrees of freedom, 
[, ]—>—*[, ] c i- In the model Eq. (2) spin operators Sj can 
be associated with classical vectors Cj and the boson op- 
erator b with a classical field a. By analogy, the Poisson 
brackets between the classical variables correspond to an 
angular momentum [C a , C^] , = — e Q ^ 7 C 7 and a boson 
field [ a, o*] c ; — i commutation relations. The Hamilton 
equations of motion for such a classical model are equiv- 
alent to the mean-field approximation. In this way Eqs. 
(24, 25) can be interpreted as the classical limit of Eqs. 
(22, 23). 

In Sec. II we wrote down explicit expressions for 
the quantum dynamics of the initial state |M;4- • • • i)- 
The initial condition for the dynamical mean-field equa- 
tions, which corresponds to this state, is all spins down, 
Cj (0) = (0,0,-1/2), and a finite amplitude of the 
bosonic field, a (0) = VaJ. Below, we solve set of differ- 
ential Eqs. (24, 25) for the evolution of classical variables 
with this initial condition. 

For a small number of excitations in a system with 
many spins, M <C N, these equations were solved in [20]. 
Using the approximation C| (t) w —1/2, the equation for 
(t) drops out from Eq. (24) and the remaining system 
of equation is harmonic. In the regime 3> A this gives 
the following solution for the bosonic mode 



a(t) 



'M cos 



(yVNtj 



(26) 



The period of the oscillatory function in this limit is 
T = 2%/(vVN^, see Figure 5 for M = 10. In the 

following will use the Rabi frequency 57 = VyN which 
was introduced in Sec. II. 

For arbitrary M the solution to Eq. (24, 25) are hy- 
perelliptic functions 35 . Here we will not consider the full 
analytic form of general solution but will use only the 
spectral analysis developed in Ref. [36]. Introducing the 
following vector function (Lax vector) of an auxiliary pa- 
rameter it, 



L(«) 



asM 

V 

V 
u 

'W 1 



E 



Cj (0) 



(27) 
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FIG. 1. Numerical solution of the equation \J L 2 {u) — with 
L(u) from Eq. (27) for N = 50, M = N, and ft/A = 10. The 
roots that merge into a continuous line are plotted in green 
and two pairs of discreet roots are plotted in blue. 



the frequency spectrum is related to the roots of the equa- 
tion yl L 2 (u) = 0. We analyze them numerically in the 
limit N 3> 1, and find that, in the high coupling regime 
fl 3> A, all roots merge into a continuous line except two 
complex conjugated pairs, see example in Figure 1. Note 
that all coefficients of the polynomial 1? (it) are real thus 
every complex root has a complex conjugated partner. 
The dynamical variables that correspond to the continu- 
ous band form a decay part of the solution and the two 
discreet frequencies give an oscillating part that we will 
be interested in. The discreet roots can be found by turn- 
ing the summation over j in Eq. (27) into an integral, 

J2j ~^ s Ju-a/2 Then, the equation sj L 2 (u) = 
turns into 



±2iga (0) 



V 2 N 
A 



In 



A/2 



u- A/2 



(28) 



where the total width of splittings A = ejv — t\ is finite. 

In the limit M = N, opposite to M <C N, the roots of 
Eq. (28) have zero real part. Parametrizing the roots as 
it = iuqA/2 we obtain 



±Wa (0) 



2V 2 N 



u 



(ir — 2 tan uq) 



(29) 



Then a 1/iio-expansion gives the imaginary parts as 
Ul , 2 = u A/2 = ±(^±^ 3 

In the case of two discreet roots the hyperelliptic func- 
tion of many variables reduces to an elliptic function 
of only one variable 36 which corresponds to an effective 
model of a single collective spin coupled to a boson. Fol- 
lowing the procedure of constructing a 1-spin solution in 
Ref. [35] we write 



(30) 



6 



iau 



(31) 



where Q4 (u) — (u 2 + uf J (u 2 + u^j is a polynomial given 
by the imaginary roots of Eq. (28). Here we choose 
u 2 > u x . 

The differential Eq. (30) defines a Jacobi elliptic func- 
tion 



u (t) = iuisn Qu^l t — A,k) , 



(32) 



where k = \u\/u2\ is the elliptic modulus, and A is an 
unknown constant of integration. Integrating the second 
equation separately for a and u we get 



a(t) = B 



dn (A) - Vkcn (A) 



dn [u<xt — A) — spkcxx (u$t — A) 



(33) 



where B is a second constant of integration. 

From the initial conditions, the phase of the oscillation 
at t = is A = 0. The second constant of integration 
is obtained from the condition a (t = 0) = Va" as B = 
y/~N. Finally, expanding the parameters in Eq. (33) for 
A 2 / (V 2 N) < 1 we obtain 



a(t) 



2A 
V3V 



(ill 



\fkcn. ^ 



V^Nt 



(34) 



The period of the oscillatory function for this initial con- 
dition is given by the complete elliptic integral of the first 

kind, T = 8K (k) I (V\fif) , see Fig. 5 for M = 50. 
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FIG. 2. Time evolution of bosonic occupation number, for 
N — 50 spins and M = 1,25,50 excitations, Coupling 
strength is fi/A = 0.002. The truncation error is 8 = 0.05%, 
using only states with a single spin excitation. 



IV. RESULTS 



Following the different regimes of interaction strength, 
we discuss now the effect of a large number of excita- 
tions on the time-evolved bosonic occupation. When the 
interaction strength is very weak, the discreetness of the 
spin subsystem plays a major role. Indeed, when the 
Rabi frequency is smaller than the spacing between 
e^s, only spins which are very close to resonance with the 
bosonic mode are significantly hybridized with the cav- 
ity while the rest only plays a weak perturbative role. 
The resulting dynamics are therefore expected to exhibit 
non-universal behavior linked to the specific choice of the 
band's tj. This results only in weak oscillations around 
(b^b^ = M. A large M (25 and 50 are shown) does not 
bring any major qualitative changes to the time evolu- 
tion, see Fig. 2. 

Considering that the non-universal behavior of this 
particular regime is characterized by only a few effective 
degrees of freedom regardless of the number of excita- 
tions, we move away to Rabi frequencies larger than the 
level spacing but still significantly smaller than the total 
bandwidth, (A/N) < O < A. Since more and more spins 
get significantly mixed with the bosonic mode, this ulti- 
mately leads to the "universal weak-field regime" which 



exhibits a decay of the bosonic population. However, 
in this particular regime, any truncation of the Hilbert 
space leads to an important loss of information and there- 
fore to a large error evidenced by a badly saturated sum 
rule (21). In trying to address the validity of the mean- 
field approach at such Rabi frequencies we therefore have 
to limit ourselves to small system sizes. In fact, the num- 
ber of eigenstates has to be small enough to be able to 
compute every one of them in a reasonable amount of 
time therefore accessing exact results. Fig. 3 presents 
a comparison of the mean-field and the exact quantum 
dynamics for a small system containing only A^ = 12 spin 
degrees of freedom. 

In the presence of a single excitation M = 1, the 
bosonic occupation number (b*b)(t) rapidly decays to al- 
most 0, which is remarkably well captured by the solu- 
tion of the semi-classical Eqs. (23). The origin of this 
decay lies in the significant overlap of many eigenstates 
of the system with the initial state, in contrast to the lim- 
ited number of important eigenstates in the ft ^> ( A /N) 
regime. Decomposing the dynamics into a persistent os- 
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FIG. 3. Time evolution of bosonic occupation number, exact 
(red solid) and meanfield (black dotted) for N = 12 spins 
and M = 1, 9 excitations. For M = 1, the decay part of eq. 
(35) is shown (green dotted-dashed). Coupling strength is 
Q/A = 0.3. 



cillation [Eq. (26] and a decay part, 



(4n 2 /A 2 ) cos (y At/ '2) dy 

(«-^Ks)) 2+ (^) s 



(35) 



as in Rcf. 19, wc indeed find that the the continuum 
part of the energy spectrum is dominant. Moving away 
from M < iV we look at a strongly excited initial state 
containing M = 9 bosons. Once again the initial decay 
is perfectly reproduced by the mean-held treatment in 
spite of the relatively large number of excitations. At 
later times the exact quantum treatment deviates from 
the mean-held approach due to small size of the system. 

Within the restriction to short times imposed by fi- 
nite size effects, the mean-held approach remains valid 
even at large excitation numbers in this regime. How- 
ever, the severe limitation on the system size makes it 
difficult to reach a precise conclusion for any larger sys- 
tems. We therefore turn our focus to the strong cou- 
pling regime, > A. In this limit, a drastic trunca- 
tion becomes possible while maintaining a satisfying sat- 
uration of the sum rule (21). We first present in Fig. 
4, the energy spectrum characterizing the time evolu- 
tion of the bosonic occupation numbers for a system of 
N = 50 spins. Specifically, we plot the work distribution 
P(E) = Eti K M ;+ \&)\ 2 5{E-Ei) which, accord- 
ing to eq. (13), describes the frequency content of the 
initial condition. For these calculations the dimension 
of the Hilbert space is reduced from d = O(10 15 ) (for 
M = N) to d = N + 1 by keeping only the states with M 
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4. Work distribution function for N = 50 spins and 
10,40,48,50 exctitations. Coupling strength is Q/A = 



divergent rapidities. Doing so maintains the truncation 
error in (Wb)(t)/M below a maximum of 5 = 5%. 

For low excitation numbers, the spectrum presents it- 
self as a series of nearly equally spaced peaks. According 
to Eq. (13), this constant energy difference between the 
contributing eigenstates indicates a periodic oscillation 
in the time evolution. However, as the number of excita- 
tions is increased, deviations from the harmonic progres- 
sion become more and more important and is particu- 
larly evident in the M = N results where the low energy 
contributions are much closer than the high energy ones. 
Additionally, the shape of the distribution is severely al- 
tered. Therefore, while one can expect the periodicity 
of the semi-classical results to be adequately reproduced 
for small M <C N, the opposite regime will be character- 
ized by a set of incommensurate frequencies ultimately 
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leading to some decay. 

This is evidenced by looking at the explicit time evolu- 
tion of the average bosonic occupation presented in Fig. 
5 for a variety of initial number of excitations. Both the 
mean-field behavior (black) and the quantum evolution 
(red) are plotted for TV = 50 spins. 




At/2 



FIG. 5. Time evolution of bosonic occupation number, exact 
(red solid) and mean-field (black dotted) for N — 50 spins 
M = 10, 40, 48, 50 and fi/A = 10. 

Two main differences between the quantum and mean- 
field dynamics are seen. First, the quantum oscillation 
frequency is systematically shifted to higher frequencies 
and, at the same time, the amplitude is shown to de- 
cay. When the number of excitations is small enough the 
frequency shift is small and the decay is slow compared 
to the time scale set by the oscillation frequency. How- 
ever, when M and N become comparable, the quantum 
dynamics exhibits a rapid amplitude decay which, not be- 
ing captured by the mean-field analysis, leads to a dras- 



tic difference between both descriptions of the bosonic 
occupation. As evidenced by Fig. 5, a larger difference 
in the oscillation frequency also occurs, making the dis- 
tinction between both approaches even more important. 
Nonetheless, even for M = N the mean-field approach 
is shown to capture perfectly the initial instability and 
provides an accurate description up to some finite time. 

In order to characterize the regime of validity of the 
mean-field approximation, we extract an Ehrenfest time 
by looking for the earliest time where the mean-field and 
quantum bosonic occupations differ by 10% of the initial 
population. The resulting times are plotted as a function 
of the excitation number in Fig. 6. 
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FIG. 6. Outer plot: Ehrenfest time for N=50 as a function 
of excitation number M for high filling factors. Dashed line 
is a t e (M) = (too + be~ cM ) fit. Inset: The same plot for the 
whole range of M. 

For few excitations M <C N both the mean-field and 
quantum numerical calculations were shown to coincide 
up to l/N corrections 19 . Here we see that the Ehren- 
fest time initially undergoes a rapid, seemingly linear, 
decrease as M increases. When the strongly excited 
regime is reached, at M « 0.87V, this behavior is dras- 
tically modified. It is then well described by an ex- 
ponential fit t e (M) — (too + be~ cM ) with parameters 
= 0.15,6 = 1003.65, c = 0.28. The saturating de- 
crease of the Ehrenfest time indicates that the mean-field 
description retains its validity in the description of the 
initial stages even as M reaches TV. 

For a low number of initial excitations 19 , the mean- 
field approximation is know to be exact in the limit 
TV — >• oo. Moreover, as evidenced in Ref. [19] and in this 
work, a modest mesoscopic number of spins TV w O(10 2 ) 
is sufficient for the mean-field treatment to adequately 
describe the dynamics over many oscillation periods. For 
a strongly excited system, we plot in Fig. 7 the system 
size dependence of the Ehrenfest time obtained at large 
fillings (M = TV, M = 0.97V). Even for a large number of 
excitations, the growth of tg with increasing system size 
indicates that the mean-field approach could still provide 
an adequate description, even at long times, for thermo- 
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V. CONCLUSIONS 

By exploiting the quantum integrability of the Dicke 
model, we were able to calculate the quantum dynamics 
of an initially populated single bosonic mode interacting 
with an ensemble of inhomogeneous ensemble of two- level 
systems. For strong enough couplings, this method based 
on the algebraic Bethe ansatz provides a simple trunca- 
tion scheme, which allowed us to treat relatively large 
systems even for a strongly excited initial state. 

We compared the numerical solutions of its non- 
equilibrium dynamics with its mean-field description. Fo- 
cusing on the strong coupling regime, where mean- field 
theory predicts oscillating periodic solutions, we confirm 
that at low excitation numbers both solutions agree up to 
a relatively long finite Ehrenfest time. However, going to 
more strongly excited systems leads to a rapid shortening 
of the mean-field description's period of validity due to a 
shift in the oscillation frequency combined with a decay 
of the oscillation's amplitude which are exclusively cap- 
tured by a full quantum treatment . For relatively large 
mesoscopic systems, we demonstrate that, although it 
cannot capture the long time dynamics, the initial decay 
of the bosonic excitations is still adequately described by 
the classical mean-field theory. 



dynamically large N —¥ oo systems. However, since this 
growth is extremely slow and therefore in stark contrast 
to the M <C N case, mesoscopic systems remain too 
small for the classical mean-field treatment to describe 
correctly the behavior past the initial stage of decay. 
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